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Abstract 



1 Introduction 



Point pattern sets arise in many different areas of 
physical, biological, and applied research, represent- 
ing many random realizations of underlying pattern 
formation mechanisms. These pattern sets can be 
heterogeneous with respect to underlying spatial pro- 
cesses, which may not be visually distiguishable. This 
heterogeneity can be elucidated by looking at statisti- 
cal measures of the patterns sets and using these mea- 
sures to divide the pattern set into distinct groups 
representing like spatial processes. We introduce here 
a numerical procedure for sorting point pattern sets 
into spatially homogenous groups using Functional 
Principal Component Analysis (FPCA) applied to 
the approximated Minkowski functionals of each pat- 
tern. We demonstrate that this procedure correctly 
sorts pattern sets into similar groups both when the 
patterns are drawn from similar processes and when 
the 2nd-order characteristics of the pattern are iden- 
tical. We highlight this routine for distinguishing 
the molecular patterning of fluorescently labeled cell 
membrane proteins, a subject of much interest in 
studies investigating complex spatial signaling pat- 
terns involved in the human immune response. 



1.1 Motivation: Why study point pat- 
terns? 

Spatial points patterns naturally arise in many ar- 
eas of research in both the physical and life sciences, 
including ecology [TJ [5], crime statistics [3j, epide- 
meology [4], economics [5], seismology [6], material 
science [7], and astronomy [5]. Whether the points 
represent molecules, trees, cell phone users, or entire 
galaxies, the spatial distributions of point patterns 
belie the underlying stochastic processes that govern 
the pattern's formation. 

A new area of point pattern analysis involves 
studying the molecular patterning of proteins on the 
surfaces of cells. Due to photo-activated localiza- 
tion microscopy (PALM) [5], a new super-resolution 
microscopy technique, cell biologists are now able 
to measure the spatial distribution of fluorescently- 
tagged membrane proteins and determine the re- 
sponse of the molecules to different stimuli (see fig- 
ure [1} . By fixing the cells on a slide and exposing 
them to laser light, researchers can activate molecules 
one by one in multiple cells, locating the center to 
within 20 nm. This new technique has resulted 
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in a wealth of new point pattern data represent- 
ing different molecules and surface treatments, and 
quantitative analysis of these patterns can contribute 
much to understanding protein-protein and protein- 
membrane interactions |1Q[ 111] . From a theoreti- 
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Figure 1: Using photoactived localization mi- 
croscopy, the fluorescenctly labeled proteins are local- 
ized by fitting a point spread function to the stochas- 
tically photoactived molecules; the final pattern rep- 
resents thousands of fhioresecent images 

cal standpoint, each pattern is a pure realization of 
an underlying spatial process and can be used to 
characterize that process. From a practical perspec- 
tive, however, it takes many experimental realiza- 
tions with finite systems to discern the underlying 
structure. Furthermore, if the point interactions are 
complex or the patterns are formed in complicated 
environments (such as the membrane of a cell), the 
amount of data needed to confidently quantify a pro- 
cess becomes large and cumbersome to analyze. This 
gives rise to the need to be able to confidentely di- 
vide large sets of patterns, sorting the patterns into 
smaller, homogenous groups that can be analyzed 
further. In addition to simplifying analysis, this type 



of sorting can also provide researchers with quick in- 
formation about the homogeneity of a process and 
the experimental parameters affect this homogeneity. 

1.2 Current Methods for Sorting Pat- 
terns 

The standard method of sorting pattern sets is as 
follows: For each pattern, one calculates a list of nu- 
merical summary characteristics (e.g. Index of Dis- 
persion, Clark-Evans index). These can be regarded 
as the "coordinates" of a pattern, to which distance- 
based clustering algorithms can be applied [12J. This 
approach presents the researcher with the task of de- 
ciding which characteristics to use, how to compare 
them (normalizing, z-cores, etc) as well as how many: 
too few may result in missing information, too many 
could result in redundancy. This adds nuance to the 
sorting, limiting the statistical conclusions that can 
be drawn, and making trustable automation of the 
sorting procedure for large pattern sets difficult to 
accomplish. 

A more robust sorting technique has been devel- 
oped, where patterns are sorted using Functional 
Principal Component Analysis (FPCA) on smoothed 
2nd-order functionals of the patterns. This routine 
treats the point set as a set of functions, {a,i(t)}, such 
as the pairwise correlation function, g(r) [13] . The 
coordinates of each pattern are then calculated by 
finding the eigenfunctions and corresponding eigen- 
values of the equation: 



v(s,t)wi(s)ds = XiWi(t). 



(1) 



Here v(s, t) is the variance-covariance function of the 
set of functionals a(r), defined as 



v{s,t) = (7V-l)- 1 ^(a J -a l (i))(a l ( S )-a l ( S )). (2) 
i=0 

The "score" of the z-th pattern on the j-th princi- 
pal component is then J a,i(t)wj(t)dt (see |14l 115]). 
Like standard PCA, the eigenvalues form a positive 
decreasing set whose truncated sum represent the to- 
tal variance encapsulated in the included principal 



2 



components. For automation, one can simply set a 
threshold for the amount of variance to be included, 
which in turn prescribes the number of coordinates 
to be used. This feature removes the arbitrariness 
of sorting patterns via the standard method, making 
FPCA an easily automatable way of quantifying the 
difference between patterns. 

However, spatial processes can create patterns with 
more structure than 2nd-order functionals can mea- 
sure. The Neyman-Scott process (NS), introduced 
to study galaxy clustering, involves randomly dis- 
tributed parent points generating clusters of varying 
size. The complexity of the parent/daughter inter- 
action gives rise to families of NS processes with the 
same pairwise correlation function [THl El [T^], de- 
spite underlying spatial differences in the patterns. 

Baddeley and Silverman also introduced a cell pro- 
cess which is built by partitioning a domain into cells 
of equal size which are then filled with a varying num- 
ber of uniformly distributed points. Though the pro- 
cess is rather regular, they showed analytically that 
their process was indistinguishable from a Poisson 
process when only considering 2nd-order functionals 
of the pattern [18], meaning that higher-order func- 
tionals must be used to resolve this ambiguity. 

1.3 Prom points to discs 

In this paper, we apply the proximity measure of 
FPCA to the approximated Minkowski functionals 
of point patterns [19]. These functionals are calcu- 
lated by centering a disc on each point and analyz- 
ing the topology of this secondary pattern of over- 
lapping discs as a function of the radius. Since the 
overlap can be very complex, involving all possible 
combinations of individual points, these functionals 
depend on all orders of interaction simultaneously. 
This makes them a more complete "fingerprint" for 
pattern comparison [HE [5D] . These functionals have 
enjoyed marked success in astrophysics [21] . soft mat- 
ter [22] . and fluid turbulence |23) . 

For completeness, we first explain the Minkowski 
functionals and how they are applied to point pattern 
analysis. We then demonstrate the sorting procedure 
by clustering sets of patterns of both synthetically 
generated data as well as biological data represent- 



ing the spatial distributions of membrane proteins. 
Using both agglomerative and divisive clustering al- 
gorithms, we show that this procedure outperforms 
FPCA clustering with 2nd-order functionals, and in 
general we demonstrate it to be a viable method for 
automatically sorting point pattern sets. 

2 Outlining the Procedure 

2.1 Minkowski functional analysis of 
point patterns 

The first step in 2-D Minkowski functional analysis 
[T§1 |2"U] is to turn a point pattern into a "secondary 
pattern" by centering a disc of radius r at the center 
of every point (see figure [5]). 

If the radius is large enough, some of these discs 
will overlap. By combining the overlapping discs, a 
pattern of differently shaped objects is formed. The 
total area, A, of this collection of objects is then just 
the total area of the discs excluding any overlapping 
area. This is the first Minkowski measure. The sec- 
ond Minkowski measure, the total perimeter, P, of 
the pattern is the perimeter of all of the shapes, which 
is again reduced from the perimeter of the individual 
discs because of overlaps. The Euler number, x, is the 
final Minkowski measure, defined as the total num- 
ber of distinct shapes or components in the window 
(created by the overlapping discs) minus the number 
of holes. 

By calculating each of these measures first for small 
radii, where the discs do not overlap, and growing the 
radius after each calculation until the entire pattern 
window is covered, the three Minkowski functionals 
A(r),P(r),x{ r ) are approximated. Because at each 
radius, the Minkowski measures depend on the loca- 
tions of all of the points simultaneously, these func- 
tionals include information about every type of spa- 
tial structure present in the pattern, completely char- 
acterizing it (a consequence of Hadwiger's Theorem 
from integral geometry, see [M|)- This feature makes 
the Minkowski functionals a more complete measure 

1 In this paper, we will deal only with two dimensional pat- 
terns, but our procedure is easily generalizable to patterns of 
any dimension. 
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Figure 2: (a) The Minkowski functionals are calcu- 
lated by imposing discs on the point pattern. This 
new secondary structure can be characterized using 
topological measures, which vary for different radii 
(b) The three reduced Minkowski functionals for a 
2-D Poisson (random) process. These functionals are 
unitless due to the normalization by the same mea- 
sure one would expect for a set of non-overlapping 
discs 



of the underlying point interactions, including infor- 
mation from all possible groupings of points. 

When comparing patterns, one actually uses the re- 
duced Minkowski functionals, namely the Minkowski 
functionals for the pattern divided by what is ex- 
pected for a set of non-overlapping discs. These are 
given by 



a(r) 
p(r) 
e(r) 



nNr 2 
P(r) 
2irNr 

N 



(3) 
(4) 
(5) 



The functionals for a Poisson process are shown in 
figure [5]b. The analysis in this paper relies exclu- 
sively on these reduced functionals, so we will not 
differentiate between the two. 

2.2 Sorting the patterns 

Our aim is to automatically sort patterns by perform- 
ing FPCA on their approximated Minkowski func- 
tionals, clustering the patterns with their individual 
scores on the principal components. We will do the 
same with the pairwise correllation function so that 
we can directly compare our method with that of 
|13j . For each pattern set, we will use enough prin- 
cipal components to account for 95% of the varia- 
tion. For the Minkowski functionals, we will calcu- 
late the principal component scores individually for 
the area, perimeter, and Euler number and then con- 
catenate the scores into a larger vector. Then, we will 
use these scores as coordinates, applying two different 
clustering algorithms: 

• Ward's method [25]: An agglomerative technique 
which seeks to minimize the total intercluster 
variance of the distances between objects. We 
chose this method because it is well known to 
the pattern analysis community, and allows us 
to directly compare our method with that of II- 
lian et al |13) . 

• Fast Weighted Modularity [55l[57]: To implement 
this routine, we first calculate the pair-wise Eu- 
clidean distance between all patterns, Dij, and 
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transform our pattern set into a weighted graph 
with edge-weights 



Wij = max(Dij) — Di 



(6) 



Then, this algorithm aims to maximize "modu- 
larity" of this weighted network, by dividing the 
set into groups where the average edge weight 
between members of the same group is higher 
than between members of different groups. We 
chose this method because it is a global cluster- 
ing routine with large popularity in the cluster 
analysis literature, and because the software im- 
plementation is able to work with very large data 
sets (millions of objects). 

By utilizing both cluster analysis algorithms, we 
can verify whatever results we obtain, and more 
completely demonstrate the efficacy of our sorting 
method using the Minkowski functionals. 



3 Testing our sorting method 

We now apply this procedure to three different data 
sets. These sets of patterns highlight three possi- 
ble situations in which one would sort patterns: 1) 
Comparing different systems, 2) Varying a parame- 
ter in an experiment, and 3) comparing different com- 
ponents of a bi-disperse system, (see Supplemental 
Material below for specifics regarding software and 
computational methods used) . Each set is comprised 
of two groups which have an a priori cluster struc- 
ture. We then apply the sorting technique, which 
allows us to calculate the percentage that is misclas- 
sified (%MC) by looking at the fraction of patterns 
that are assigned to a group that is dominated by a 
different pattern type. 



3.1 Data Set 1: 

cesses 



Two Strauss Pro- 



The Strauss process [28\ is a germ-grain pattern sim- 
ulation model specified by two parameters, a radius 
r £ M. + and an interaction parameter 7 £ [0, 1]. The 
interaction parameter determines if grains of radius 



r will be allowed to overlap during the formation of 
the pattern (see figure [3] a). 

If 7 is small, there is a strong repulsion between 
grains, where 7 = yield a hard-core process. If 7 is 
close to 1, the repulsion is weak, where 7 = 1 yields a 
completely random process. In [13] . it was reported 
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Figure 3: (a) The regularity of a Strauss process is 
completely determined by 7, the interaction parame- 
ter (b) The left most panes display the g(r) and x( r ) 
for the 40 simulated Strauss processes. On the right 
are the results of using FPCA scores to divide the pat- 
tern set into two groups. As can be seen, both g(r) 
and the Minkowski functionals can perfectly separate 
the set into two groups corresponding to different val- 
ues of 7. 

that even for comparison of pattern sets with sim- 
ilarly strong repulsion (7 = 0.0 and 7 = 0.1) the 
pairwise correllation function was able to effectively 
distinguish different Strauss processes. We here re- 
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peat this test with 20 patterns each, fixing the num- 
ber of points at N = 1000, and letting r = 0.025. We 
also fix the number intensity to be A = (2r) -2 , which 
forces interaction between the points. 

As can be seen in figure [3]b, both g{r) and the 
Minkowski functionals are able to distinguish the two 
Strauss processes, separating the pattern set into two 
homogenous groups. This is to be expected for g(r), 
as 2nd-order interactions dominate the process, and 
is consistent with the findings of |13] . 

3.2 Data Set 2: Baddelley-Silverman 
vs. Random 

The Baddelley-Silverman process |18| is built by par- 
titioning a domain into a grid and moving from box 
to box, distributing N points in each box uniformly. 
N itself is a random number, taking on the values 0,1, 
and 10 with probabilities 1/10, 8/9, and 1/90, respec- 
tively. This causes the process to be rather regular, 
but with some strong clustering occuring every now 
and then (see figure 6). 

Since E[N] — Var[A] = 1, it can be shown that the 
Baddelley-Silverman process shares all of the same 
2nd-order characteristics as a random (CSR) process. 
In ref . [13] , the Minkowski functionals were shown to 
be able to distinguish these two processes. Therefore, 
we expect to see proper sorting when using a(r), p(r), 
and e(r), and failure using g(r). 

We simulated 29 BS processes and 29 Poisson pro- 
cesses, fixing the point number N = 1024. Using both 
the pairwise correlation function and the Minkowski 
functionals, and sorted the pattern sets to into two 
groups. As can be seen in figure [4]b, the pair- 
wise correllation function fails to sort the patterns 
correctly, creating two heterogenous groups (%MC 
> 40). However, the Minkowski functionals suc- 
cessfully divide the pattern set into two homogenous 
groups. 

3.3 Data Set 3: Bi-disperse patterns 
of inter-cellular proteins 

For an application to an experimental data set, we 
look at super-resolution images of two proteins re- 
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Figure 4: (a) A Baddelley-Silverman process side-by- 
side with a Poisson process. Despite the visible differ- 
ences, the pairwise correlation functions are identical 
(b) The left most panes display the g(r) and x(r) 
for the 58 patterns simulated. On the right are the 
results of using FPCA scores to divide the pattern 
set into two groups. As can be seen, FPCA sorting 
with g(r) creates two perfectly heterogenous groups, 
while FPCA sorting with the Minkowski functionals 
groups the patterns correctly 



G 



siding at the membrane of immune cells (see figure 
la). 

One protein under study is LAT, short for "Linker 
for Activation of T-cells", a naturally occuring pro- 
tein crucially involved in the reactions that regulate 
T-Cell antigen-dependent activation, a critical event 
in the adaptive immune response. LAT proteins have 
been seen to form clusters on the membrane with po- 
tentially complicated hierarchies [TT]. However, the 
membrane of the cell can have a lst-order effect on 
the molecular patterning of membrane proteins. It 
has been found that the location of other membrane 
protein clusters often correlates with how close the 
membrane is to the surface, and anti-correlates with 
regions of high membrane fluctuations [35] . 

Another protein, TAC (the alpha chain of the IL- 
2 receptor), can also be localized and differentiated 
from LAT by tagging with a different fluorescent 
molecule and using two different lasers with differ- 
ent wavelengths. TAC is a membrane protein that 
does not form clusters, instead distributing uniformly 
in regions where protein-membrane interactions have 
not excluded proteins. This means that TAC can 
serve as a membrane marker when studying the clus- 
tering of other proteins. Since LAT and TAC are part 
of separate signalling pathways, they also do not in- 
teract biochemically |30j . Therefore, upon sorting, 
we should get two homogenous groups representing 
the two different molecules. 

Applying FPCA on the approximated pairwise cor- 
relation functions of these data sets again yields 
strongly heterogeneous groups (%MC « 50%). This 
is visible in the pair-correlation functions (figure[5]b), 
where the individual patterns exhibit large variabil- 
ity. In contrast, because the Minkowski function- 
als consider more than just 2nd-order interactions, 
the Euler number is able to visibly distinguish the 
molecular patterns, and the pattern sorting is im- 
proved (%MC « 25%). Further success is achieved if 
we look at how FPCA sorting with each functional 
performs on its own (see figure [5]). When only us- 
ing the area, the sorting is identical to the sorting 
based on all three Minkowski functionals. However, 
sorting the LAT/TAC protein pattern set improves 
when just using the perimeter, and we achieve almost 
perfect classification when using the Euler number, 
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Figure 5: (a) The two proteins are both dispersed in 
the cell membrane, but can be visualized separately 
(b) The left most panes display the g(r) and x( r ) 
for the 16 molecular patterns. On the right are the 
results of using FPCA scores to divide the pattern 
set into two groups. As can be seen, using Minkowski 
functionals with FPCA improves the differentiation 
of the two sets 
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Figure 6: Considering the individual functionals 
(sorting using weighted modularity), the Euler num- 
ber outperforms the other two, only misclassifying 
one pattern 
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only misclassifying one pattern. This is not surpris- 
ing, since the area and perimeter are constrained 
to be smooth, positive, and monotonically decreas- 
ing, and thus cannot vary as much while the Euler 
number can vary more wildly. 



4 Conclusions and Discussion 

In this work, we have introduced the procedure to au- 
tomatically sort point pattern sets using the approxi- 
mated Minkowski functionals and Functional Princi- 
pal Component Analysis (FPCA). Using Strauss pro- 
cesses with strong repulsion, we have shown that this 
method can accurately sort point pattern sets drawn 
from very similar processes. Further, this method 
also distinguishes Baddeley-Silverman processes from 
Poisson processes, a task which the pairwise correla- 
tion function perfectly fails to accomplish. We then 
found that when looking at experimental point pat- 
terns representing proteins, FPCA sorting using the 
Minkowski functionals outperformed FPCA sorting 
with the pairwise correlation function. This suffi- 
ciently demonstrates that the Minkowski functionals 
can successfully quantify the differences between pat- 
tern sets showing complex behavior. 

We also found that FPCA sorting using only the 
Euler number strongly outperforms the other two. 
While mathematically the three functionals do com- 
pletely classify a pattern, the area and perimeter 
may only be slightly different for different spatial 
processes. This means that error introduced when 
approximating the functionals numerically may blur 
these differences, resulting in improper sorting. Since 
the Euler number is allowed to vary more dramati- 
cally as the discs combine and holes form, it can visu- 
ally distinguish very similar pattern sets, and there- 
fore leads to better sorting. 

Though we have presented this technique as a way 
to sort patterns into distinct sets for further analysis, 
the sorting itself can serve as an analysis tool. We are 
currently working to apply this tool to examine how 
the presence of different chemical cues effect the clus- 
tering of LAT proteins, as well as how T-cell activa- 
tion perturbs the patterns. Because of the Minkowski 
functional's ability to robustly characterize a pattern, 



we can treat the membership of a pattern in a partic- 
ular group as a sign of similarity between it and it is 
co-members. We can therefore look at group statis- 
tics to determine what experimental variables change 
the molecular patterns, and to what degree, allowing 
for systematic large-scale investigations of the mem- 
brane proteins and their response to different stimuli. 
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6 Supplemental Material 
6.1 Software 

To approximate the 2-D Minkowski functionals 
of our patterns, we relied exclusively on the soft- 
ware described in [19J, which was available online at 
http : //www.mathe . tu-freiberg.de/ inst/ stoch/Stoyan/morph2D, 
El This program takes as input rmin, dr, and 
rmax. Since our interest is in automation we 
used the same values for all of our patterns 
(r min =dr = M,r max = 100). 

For both smoothing and applying the Functional 
Principal Component Analysis, we used the Func- 
tional Data Analysis MATLAB packages that are 
available online at www.functionaldata.org, and 
we relied on their description in [15 for implemen- 
tation. Mimicking the procedure of [13], we first 
smoothed our functionals using cubic b-splines. 

To cluster using Ward's method, we first uti- 
lized MATLAB 's implementation in their "linkage" 
function. To implement modularity maximiza- 
tion, we used the weighted version of the Fast 
Modularity algorithm which can be found online at 



1 At the time of this paper's submission, this website was 
down; we are in the process of notifying the appropriate people 
about this issue 
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http://cs .unm.edu/~aaron/research/f astmodularift.j4htnl :> attern Simulation 



The specifics of the algorithm are the same as in [26] , 
but maximizes the weighted definition of modularity 
(for a description of this alteration, see [27]). 

Other home made programs were written to com- 
pute the second order functionals, simulate point pat- 
terns, and implement various portions of the project 
(either in MATLAB or C). Those interested in dis- 
cussing these programs should contact the authors. 



6.2 Intensity Scaling 

As reported in [24|, the Minkowski functionals are 
homogenous with regard to domain scaling. To be 
specific, for any parameter A > and domain fl C 
M d , the n.-th Minkowski functional M n (|fi|) satisfies 
the relation 



M„(|Afi|) = x d - n M n (\n\) 



(7) 



This means two patterns with different overall num- 
ber intensity will have different Minkowski function- 
als even if they are the same type of pattern. To 
address this in our pattern comparison, all patterns 
were scaled to unit intensity before their Minkowski 
functionals were approximated. 



6.3 Approximating g(r) 

To approximate the pairwise correlation function 
g(r), we used the estimator 



N 

g(r) = N-'X- 1 w iHr, dr) ^ - ?j \ < r) 

i=l 



(8) 

Here, I(x) is the indicator function and A is the num- 
ber intensity. The weight w; is the portion of the area 
of the disc centered on r\ with inner radius r and 
outer radius r + dr that is contained in the pattern 
window. We found that this method achieved better 
results than that of |13| . where g(r) is approximated 
by exploiting it's relation to the derivative of Ripley's 
K-function. 



Binomial processes were simulated using MATLAB's 
built in random number generator and scaling the 
results. MATLAB code to simulate Strauss pro- 
cesses can be found in [3T], and Baddelley-Silverman 
processes were simulated using homemade software 
based on the procedure described in [15] . 
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